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Abstract 

We show that the discrete complex, and numerous hypercomplex, Fourier transforms defined and used 
so far by a number of researchers can be unified into a single framework based on a matrix exponential 
version of Euler's formula e j9 = cos 8 + j sin 0, and a matrix root of — 1 isomorphic to the imaginary root 
j. The transforms thus defined can be computed using standard matrix multiplications and additions 
with no hypercomplex code, the complex or hypercomplex algebra being represented by the form of the 
matrix root of —1, so that the matrix multiplications are equivalent to multiplications in the appropriate 
algebra. We present examples from the complex, quaternion and biquaternion algebras, and from Clifford 
algebras Ci\ % \ and Cfe.o- The significance of this result is both in the theoretical unification, and also 
in the scope it affords for insight into the structure of the various transforms, since the formulation 
is such a simple generalization of the classic complex case, ft also shows that hypercomplex discrete 
Fourier transforms may be computed using standard matrix arithmetic packages without the need for 
a hypercomplex library, which is of importance in providing a reference implementation for verifying 
implementations based on hypercomplex code. 

1 Introduction 

The discrete Fourier transform is widely known and used in signal and image processing, and in many other 
fields where data is analyzed for frequency content pQ. The discrete Fourier transform in one dimension is 
classically formulated as: 

M-l 

M-l v ' 

/M=r]T ^>]cxp( 3 2n^) 

u=0 

where j is the imaginary root of —1, f[m] is real or complex valued with M samples, F[u] is complex valued, 
also with M samples, and the two scale factors S and T must multiply to 1/M. If the transforms are to be 
unitary then S must equal T also. In this paper we discuss the formulation of the transform using a matrix 
exponential form of Euler's formula in which the imaginary square root of —1 is replaced by an isomorphic 
matrix root. This formulation works for the complex DFT, but more importantly, it works for hypercomplex 
DFTs (reviewed in §[2|). The matrix exponential formulation is equivalent to all the known hypercomplex 
generalizations of the DFT known to the authors, based on quaternion, biquaternion or Clifford algebras, 
through a suitable choice of matrix root of —1, isomorphic to a root of —1 in the corresponding hypercomplex 
algebra. All associative hypercomplex algebras (and indeed the complex algebra) are known to be isomorphic 
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to matrix algebras over the reals or the complex numbers. For example, Ward [2J § 2.8] discusses isomorphism 
between the quaternions and 4x4 real or 2 x 2 complex matrices so that quaternions can be replaced by 
matrices, the rules of matrix multiplication then being equivalent to the rules of quaternion multiplication 
by virtue of the pattern of the elements of the quaternion within the matrix. Also in the quaternion case, 
Ickes [3j wrote an important paper showing how multiplication of quaternions could be accomplished using a 
matrix- vector or vector-matrix product that could accommodate reversal of the product ordering by a partial 
transposition within the matrix. This paper, more than any other, led us to the observations presented here. 

The fact that a hypercomplex DFT may be formulated using a matrix exponential may not be surprising. 
Nevertheless, to our knowledge, those who have worked on hypercomplex DFTs have not so far noted or 
exploited the observations made in this paper, which is surprising, given the ramifications discussed later. 

2 Hypercomplex transforms 

The first published descriptions of hypercomplex transforms that we are aware of date from the late 1980s, 
using quaternions. In all three known earliest formulations, the transforms were defined for two-dimensional 
signals (that is functions of two independent variables). The two earliest formulations [H §6.4.2] and [5J 
Eqn. 20] are almost equivalent (they differ only in the placing of the exponentials and the signal and the 
signs inside the exponentials) 1 : 



In a non-commutative algebra the ordering of exponentials within an integral is significant, and of course, 
two exponentials with different roots of -1 cannot be combined trivially. Therefore there are other possible 
transforms that can be defined by positioning the exponentials differently. The first transform in which the 
exponentials were placed either side of the signal function was that of Ell [6j [7] : 



This style of transform was followed by Chernov, Bulow and Sommer [H OH [10] and others since. In 1998 
the present authors described a single-sided hypercomplex transform for the first time [llj exactly as in 
([I]) except that / and F were quaternion- valued and j was replaced by a general quaternion root of — 1. 
Expressed in the same form as the transforms above, this would be: 



where /x is now an arbitrary root of -1, not necessarily a basis element of the algebra. The realisation that an 
arbitrary root of -1 could be used meant that it was possible to define a hypercomplex transform applicable 
to one dimension: 



Pei et al have studied efficient implementation of quaternion FFTs and presented a transform based on 
commutative reduced biquaternions [121 1 1 3j . Ebling and Scheuermann defined a Clifford Fourier transform 
[141, § 5.2], but their transform used the pseudoscalar (one of the basis elements of the algebra) as the square 
root of —1. 

Apart from the works by the present authors [TTJ [T51 [TH] , the idea of using a root of — 1 different to 
the basis elements of a hypercomplex algebra was not developed further until 2006, with the publication of 

In comparing the various formulations of hypercomplex transforms, we have changed the symbols used by the original 
authors in order to make the comparisons clearer. We have also made trivial changes such as the choice of basis elements used 
in the exponentials. 
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a paper setting out the roots of —1 in biquaternions (a quaternion algebra with complex numbers as the 
components of the quaternions) |17| . This work prepared the ground for a biquaternion Fourier transform 
[18| based on the present authors' one-sided quaternion transform [11J . More recently, the idea of finding 
roots of —1 in other algebras has been advanced in Clifford algebras by Hitzer and Ablamowicz [TH] with 
the express intent of using them in Clifford Fourier transforms, perhaps generalising the ideas of Ebling and 
Scheuermann |14j . Finally, in this very brief summary of prior work we mention that the idea of applying 
hypercomplex algebras in signal processing has been studied by other authors apart from those referenced 
above. For an overview see [2D] . 

In what follows we concentrate on DFTs in one dimension for simplicity, returning to the two dimensional 
case in §[7] 



3 Matrix formulation of the discrete Fourier transform 
3.1 Matrix form of Euler's formula 

The transform presented in this paper depends on a generalization of Euler's formula: cxp i9 — cos 9 + i sin 9, 
in which the imaginary root of —1 is replaced by a matrix root, that is, a matrix that squares to give a 
negated identity matrix. Even among 2x2 matrices there is an infinite number of such roots [211 pl6]. In the 
matrix generalization, the exponential must, of course, be a matrix exponential [22[ §11.3]. The following 
Lemma is not claimed to be original but we have not been able to locate any published source that we could 
cite here. Since the result is essential to Theorem [IJ we set it out in full. 

Lemma 1. Euler's formula may be generalized as follows: 

e je = Icos9 + Jsinfl 
where I is an identity matrix, and J 2 = —I. 

Proof. The result follows from the series expansions of the matrix exponential and the trigonometric func- 
tions. From the definition of the matrix exponential § 11.3]: 

k=0 

Noting that J° = I (see [23l Index Laws]), and separating the series into even and odd terms: 

I9 2 7(9 4 I9 6 
~^T + _ 4! 6!~ + '" 
J9 3 J9 5 J9 7 

^r + - 

= J cos 9 + J sin ( 



M ~ 3! ' 5! 



□ □ 

Note that matrix versions of the trigonometric functions are not needed to compute the matrix exponen- 
tial, because 9 is a scalar. In fact, if the exponential is evaluated numerically using a matrix exponential 
algorithm or function, the trigonometric functions are not even explicitly evaluated |22i §11.3]. In practice, 
given that this is a special case of the matrix exponential, (because J 2 = —J), it is likely to be numerically 
preferable to evaluate the trigonometric functions and to sum scaled versions of J and J. 

Notice that the matrix e has a structure with the cosine of 9 on the diagonal and the (scaled) sine of 
9 where there arc non-zero elements of J . 
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3.2 Matrix form of DFT 



The classic discrete Fourier transform of ([1} may be generalized to a matrix form in which the signals are 
vector- valued with N components each and the root of — 1 is replaced by an N X N matrix root J such that 
J 2 = I. In this form, subject to choosing the correct representation for the matrix root of —1, we may 
represent a wide variety of complex and hypercomplex Fourier transforms. 



Theorem 1. The following are a discrete Fourier transform pair 2 : 

M-l 

F[:, U ] = 5^exp(-J27r^)/[:,m] 

m=0 
A/-1 

f[:,m] =T E exp( J2n—)F[:,u] 



(2) 



(3) 



u=0 



where J is a N x N matrix root of — 1, / and F are N x M matrices with one sample per column, and the 
two scale factors S and T multiply to give 1/AI. 

Proof. The proof is based on substitution of the forward transform ^ into the inverse © followed by 
algebraic reduction to a result equal to the original signal /. We start by substituting (J2)) into ([3]), replacing 
m by M to keep the two indices distinct, and at the same time replacing the two scale factors by their 
product 1/M: 



M-l 



J2-n 



M-l 

E< 

M=0 



The exponential of the outer summation can be moved inside the inner, because it is constant with respect 
to the summation index A4: 



M-l M-l 



m e 



u=0 M=0 

The two exponentials have the same root of — 1, namely J, and therefore they can be combined: 

M-i M-l 

M 



1 M-l A/-1 



(m-MW 



f[-,M] 



u=0 M=0 



We now isolate out from the inner summation the case where m = Ai . In this case the exponential reduces 
to an identity matrix, and we have: 

M-l 

/[•m] = -^/[-m] 



u=0 
M-l 



-E 

M ^ 



M-l 

E 

A4=0 



-f[;M] 



The first line on the right sums to /[:, m], which is the original signal, as required. To complete the proof, we 
have to show that the second line on the right reduces to zero. Taking the second line alone, and changing 
the order of summation, wc obtain: 



M-i 

E 

M=0 



M=j£m 



M-l 

E< 

u=0 



J2tt± 



fi;M] 



2 The colon notation used here will be familiar to users of MATLAB® (an explanation may be found in 1221 § 1.1.8]). Briefly, 
f[:,m] means the m th column of the matrix /. 
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Using Lemma [T] we now write the matrix exponential as the sum of a cosine and sine term. 



M-l 

E 

M=0 



M-l 



/ (m-M)u 
J> cos 2tt- — — 

-J ^ sin 2tt 



f[-,M] 



Since both of the inner summations are sinusoids summed over an integral number of cyles, they vanish, and 
this completes the proof. □ □ 

Notice that the requirement for J 2 = —I is the only constraint on J. It is not necessary to constrain 
elements of J to be real. Note that J 2 = —I implies that J 1 = —J, hence the inverse transform is obtained 
by negating or inverting the matrix root of —1 (the two operations are equivalent). 

The matrix dimensions must be consistent according to the ordering inside the summation. As written 
above, for a complex transform represented in matrix form, / and F must have two rows and M columns. 
If the exponential were to be placed on the right, / and F would have to be transposed, with two columns 
and M rows. 

It is important to realize that ([2]) is totally different to the classical matrix formulation of the discrete 
Fourier transform, as given for example by Golub and Van Loan [22l §4.6.4]. The classic DFT given in (JXJ) 
can be formulated as a matrix equation in which a large M x M Vandermonde matrix containing n th roots 
of unity multiplies the signal / expressed as a vector of real or complex values. Instead, in ^ each matrix 
exponential multiplies one column of /, corresponding to one sample of the signal represented by / and 
the dimensions of the matrix exponential are set by the dimensionality of the algebra (2 for complex, 4 for 
quaternions etc.). In ([2]) it is the multiplication of the exponential and the signal samples, dependent on the 
algebra involved, that is expressed in matrix form, not the structure of the transform itself. 

Readers who are already familiar with hypercomplex Fourier transforms should note that the ordering of 
the exponential within the summation ^ is not related to the ordering within the hypercomplex formula- 
tion of the transform (which is significant because of non-commutative multiplication) . The hypercomplex 
ordering can be accommodated within the framework presented here by changing the representation of the 
matrix root of —1, in a non-trivial way, shown for the quaternion case by Ickes [3l Equation 10] and called 
transmutation. We have studied the generalisation of Ickes' transmutation to the case of Clifford algebras, 
and it appears that there is a more general operation. In the cases we have studied this can be described as 
negation of the off-diagonal elements of the lower-right sub- matrix, excluding the first row and column 3 . We 
believe a more general result is known in Clifford algebra but we have not been able to locate a clear state- 
ment that we could cite. We therefore leave this for later work, as a full generalisation to Clifford algebras 
of arbitrary dimension requires further work, and is more appropriate to a more mathematical paper. 



4 Examples in specific algebras 

In this section we present the information necessary for ([2]) and ([3]) to be verified numerically. In each of the 
cases below, we present an example root of —1 and a matrix representation 4 . We include in the Appendix 
a short MATLAB® function for computing the transform in ©. The same code will compute the inverse by 
negating J. This may be used to verify the results in the next section and to compare the results obtained 
with the classic complex FFT. In order to verify the quaternion or biquaternion results, the reader will need 
to install the QTFM toolbox [24], or use some other specialised software for computing with quaternions. 

3 This gives the same result as transmutation in the quaternion case. 

4 The matrix representations of roots of -1 are not unique - a transpose of the matrix, for example, is equally valid. The 
operations that leave the square of the matrix invariant probably correspond to fundamental operations in the hypercomplex 
algebra, for example negation, conjugation, reversion. 
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4.1 Complex algebra 

The 2x2 real matrix (j ~J) can be easily verified by eye to be a square root of the negated identity matrix 
( ~q _° ) , and it is easy to verify numerically that Euler's formula gives the same numerical results in the 
classic complex case and in the matrix case for an arbitrary 9. This root of —1 is based on the well-known 
isomorphism between a complex number a + jb and the matrix representation ( j Theorem 1.6] 5 . 
The structure of a matrix exponential e je using the above matrix for J is ( g ^ ) where C = cos and 
S = sine. 



4.2 Quaternion algebra 

The quaternion roots of —1 were discovered by Hamilton pp203, 209], and consist of all unit pure 
quaternions, that is quaternions of the form xi + yj + zk subject to the constraint x 2 + y 2 + z 2 = 1. A 
simple example is the quaternion fi = (i + j + fe)/v3, which can be verified by hand to be a square root of 
— 1 using the usual rules for multiplying the quaternion basis elements (i 2 = j 2 = k 2 = ijk = —1). Using 
the isomorphism with 4x4 matrices given by Ward [2j § 2.8], between the quaternion iv + xi + yj + zk and 
the matrix: 



f w 


—x 


-y 




X 


It' 




I) 


y 


z 


w 


—x 


V 


-y 


X 


w J 



we have the following matrix representation: 





(o 


-i 


-l 


-A 


1 


i 





-i 


i 




i 


l 





-i 




V 


-i 


1 





Notice the structure that is apparent in this matrix: the 2x2 blocks on the leading diagonal at the top left 
and bottom right can be recognised as roots of —1 in the complex algebra as shown in S I4.1I 

Proposition 1. Any matrix of the form: 



A) 


—x 


-y 




X 





—z 


y 


y 


z 





—x 


V- 


-y 


X 





with x 2 + y 2 + z 2 = 1 is the square root of a negated 4x4 identity matrix. Thus the matrix representations 
of the quaternion roots of —1 are all roots of the negated 4x4 identity matrix. 

Proof. The matrix is anti-symmetric, and the inner product of the i th row and i th column is — x 2 — y 2 — z 2 , 
which is —1 because of the stated constraint. Therefore the diagonal elements of the square of the matrix 
are —1. Note that the rows of the matrix have one or three negative values, whereas the columns have zero 
or two. The product of the i th row with the j th column, i ^ j, is the sum of two values of opposite sign and 
equal magnitude. Therefore all off-diagonal elements of the square of the matrix are zero. □ □ 

The structure of a matrix exponential e je using a matrix as in Proposition Q] for J is: 

/ C —xS —yS —zS\ 

xS C —zS yS 

yS zS C —xS 

\zS —yS xS C) 

5 We have used the transpose of Ward's representation for consistency with the quaternion and biquaternion representations 
in the two following sections. 
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where, as before, C = cos 6* and S = sin 8. 



4.3 Biquaternion algebra 

The biquaternion algebra [2| Chapter 3] (quaternions with complex elements) can be handled exactly as in 
the previous section, except that the 4x4 matrix representing the root of — 1 must be complex (and the signal 
matrix must have four complex elements per column). The set of square roots of —1 in the biquaternion 
algebra is given in [17]. A simple example is i+j + k + I(j — k) (where I denotes the classical complex 
root of —1, that is the biquaternion has real part i + j + k and imaginary part j — k). Again, this can be 
verified by hand to be a root of —1 and its matrix representation is: 



/ o 
1 

1+1 

V- 1 



o 

i - 1 

-1-7 



-I- I 
-1 + / 



1 



+ 
1+1 
-1 

o / 



Again, sub-blocks of the matrix have recognizable structure. The upper left and lower right diagonal 2x2 
blocks are roots of —1, while the lower left and upper right off-diagonal 2x2 blocks arc nilpotent - that is 
their square vanishes. 



4.4 Clifford algebras 

Recent work by Hitzer and Ablamowicz [19] has explored the roots of —1 in Clifford algebras Ci VA up to 
those with p + q = 4, which are 16-dimcnsional algebras 6 . The derivations of the roots of -1 for the 16- 
dimensional algebras are long and difficult. Therefore, for the moment, we confine the discussion here to 
lower-order algebras, noting that, since all Clifford algebras arc isomorphic to a matrix algebra, we can be 
assured that if roots of -1 exist, they must have a matrix representation. Using the results obtained by Hitzer 
and Ablamowicz, and by finding from first principles the layout of a real matrix isomorphic to a Clifford 
multivector in a given algebra, it has been possible to verify that the transform formulation presented in 
this paper is applicable to at least the lower order Clifford algebras. The quaternions and biquaternions are 
isomorphic to the Clifford algebras C£o,2 and Ci^,o respectively so this is not surprising. Nevertheless, it is 
an important finding, because until now quaternion and Clifford Fourier transforms were defined in different 
ways, using different terminology, and it was difficult to make comparisons between the two. Now, with 
the matrix exponential formulation, it is possible to handle quaternion and Clifford transforms (and indeed 
transforms in different Clifford algebras) within the same algebraic and/or numerical framework. 

We present examples here from two of the 4-dimensional Clifford algebras, namely C£x,i and C£ 2: q. These 
results have been verified against the CLICAL package [55] to ensure that the multiplication rules have been 
followed correctly and that the roots of —1 found by Hitzer and Ablamowicz are correct. 

Following the notation in [19j . we write a multivector in Ct\,\ as a + b\ei + b 2 e 2 + /3e±2, where e\ = 
+1, e\ = — 1, e\ 2 = +1 and eie 2 = ei2- A possible real matrix representation is as follows: 



fa 


bi 


-bi 




h 


a 


-P 


h 


b 2 


-0 


a 


h 




-b 2 


bi 


aj 



In this algebra, the constraints on the coefficients of a multivector for it to be a root of —1 are as follows: 
o = and b\ - b\ + f3 2 = -1 QU Table l] 7 . Choosing &i = £ = 1 gives b 2 = V3 and thus e\ + \/3e 2 + e 12 

6 p and q are non- negative integers such that p + q = n and n > 1. The dimension of the algebra (strictly the dimension of 
the space spanned by the basis elements of the algebra) is 2 n . 

7 We have re-arranged the constraint compared to 1191 Table 1] to make the comparison with the quaternion case easier: we 
see that the signs of the squares of the coefficients match the signs of the squared basis elements. 
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which can be verified algebraically or in CLICAL to be a root of —1. The corresponding matrix is then: 



(o 


1 


-V3 




1 





-1 


^) 




-1 





1 


V 1 


-V3 


1 


0/ 



Following the same notation in algebra C£2,0j m which e 2 = e\ = +1, e 2 2 = — 1, a possible matrix represen- 
tation is: 



/a 


61 






&i 


a 


P 


-b 2 






a 


h 




-62 


h 


a) 



The constraints on the coefficients arc a = and b\ + b\ — /3 2 = — 1, and choosing b\ = bi = 1 gives /? = \/3 
and thus ei + e2 + V / 3ei2 is a root of — 1. The corresponding matrix is then: 



/o 


1 


1 


-\/3\ 


1 





V3 


-1 


1 


-V3 





1 


VV3 


-1 


1 


0/ 



Notice that in both of these algebras the matrix representation of a root of —1 is very similar to that given 
for the quaternion case in Proposition [1] with zeros on the leading diagonal, an odd number of negative 
values in each row and an even number in each column. It is therefore simple to see that minor modifications 
to Proposition [T] would cover these algebras and the matrices presented above. 



5 An example not based on a specific 
algebra 

We show here using an arbitrary 2x2 matrix root of —1, that it is possible to define a Fourier transform 
without a specific algebra. Let an arbitrary real matrix be given as J = ( " b d ) , then by brute force expansion 
of J 2 = — / we find the original four equations reduce to but two independent equations. Picking (a, b) and 
solving for the remaining coefficients we find that any matrix of the form: 

( a h \ 

^-(l + a 2 )^ -a) 

with finite a and b, and b ^ 0, is a root of — 1. Choosing instead (a, c) we get the transpose form: 

-(l + a 2 )/c^ 

where c 7^ 0. Choosing the cross-diagonal terms (b, c) yields: 




where K = v— 1 — be and be < — 1. 

In all cases the resulting matrix has eigenvalues of A = ±i. (This is a direct consequence of the fact 
that this matrix squares to — 1.) Each form, however, has different eigenvectors. The standard matrix 
representation for the complex operator i is (Jo) w ith eigenvectors v = [l,±i]. In the matrix with (a, b) 
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parameters the eigenvectors are v = [1, —b/(a ± i)] whereas the cross-diagonal form with (b, c) parameters 
has eigenvectors v = [1, (k ± i)/c\. 

These forms suggest the interesting question: which algebra, if any, applies here 8 ; and how can the Fourier 
coefficients (the 'spectrum') be interpreted? We are not able to answer the first question in this paper. The 
'interpretation' of the spectrum is relatively simple. Consider a spectrum F containing only one non-zero 
column at index uq with value (y) and invert this spectrum using ([3]). Ignoring the scale factor, the result 
will be the signal: 

/tH (J) 

The form of the matrix exponential depends on J. In the classic complex case, as given in § 14. 1[ the matrix 
exponential, as already seen, takes the form: 

/cos 9 — sin^\ 
l^sin^ cos 9 J 

where 9 = 2ir I! jf l . This is a rotation matrix and it maps a real unit vector ( J ) to a point on a circle in the 
complex plane. It embodies the standard phasor concept associated with sinusoidal functions. Using the 
same analysis, this time using the matrix in Q above, one obtains for the matrix exponential the 'phasor' 
operator: 

/ cos 9 + k sin 9 b sin 9\ 

\ c sin 9 cos 9 — n sin 9 J 

Instead of mapping a real unit vector ( J ) to a point on a circle, this matrix maps to an ellipse. Thus, we 
see that a transform based on a matrix such as that in @ has basis functions that are projections of an 
elliptical, rather than a circular path in the complex plane, as in the classical complex Fourier transform. We 
refer the reader to a discussion on a similar point for the one-sided quaternion discrete Fourier transform in 
our own 2007 paper [16j § VI], in which we showed that the quaternion coefficients of the Fourier spectrum 
also represent elliptical paths through the space of the signal samples. 

It is possible that the matrices discussed in this section could be transformed by similarity transformations 
into matrices representing elements of a Clifford algebra 9 . Note that in the quaternion case, any root of -1 
lies on the unit sphere in 3-space, and can therefore be transformed into another root of -1 by a rotation. 
It is possible that the same applies in other algebras, the transformation needed being dependent on the 
geometry. Clearly there are interesting issues to be studied here, and further work to be done. 



6 Non-existence of transforms in algebras with odd dimension 

In this section we show that there are no real matrix roots of — 1 with odd dimension. This is not unexpected, 
since the existence of such roots would suggest the existence of a hypercomplex algebra of odd dimension. 
The significance of this result is to show that there is no discrete Fourier transform as formulated in Theorem 
[U for an algebra of dimension 3, which is of importance for the processing of signals representing physical 
3-space quantities, or the values of colour image pixels. We thus conclude that the choice of quaternion 
Fourier transforms or a Clifford Fourier transform of dimension 4 is inevitable in these cases. This is not an 
unexpected conclusion, nevertheless, in the experience of the authors, some researchers in signal and image 
processing hesitate to accept the idea of using four dimensions to handle three-dimensional samples or pixels. 
(This is despite the rather obvious parallel of needing two dimensions - complex numbers - to represent the 
Fourier coefficients of a real- valued signal or image.) 

Theorem 2. There are no N x TV matrices J with real elements such that J 2 = I for odd values of N. 

8 It is possible that there is no corresponding 'algebra' in the usual sense. Note that there are only two Clifford algebras of 
dimension 2, one of which is the algebra of complex numbers. The other has no multivector roots of -1 19 § 4] and therefore 
the roots of —1 given above cannot be a root of —1 in any Clifford algebra. 

9 We are grateful to Dr Eckhard Hitzer for pointing this out, in September 2010. 
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Proof. The determinant of a diagonal matrix is the product of its diagonal entries. Therefore | — I\ = — 1 
for odd TV. Since the product of two determinants is the determinant of the product, |«7 2 | = —1 requires 
| J| 2 = —1, which cannot be satisfied if J has real elements. □ □ 



7 Extension to two-sided DFTs 

There have been various definitions of two sided hypercomplex Fourier transforms and DFTs. We consider 
here only one case to demonstrate that the approach presented in this paper is applicable to two-sided as 
well as one-sided transforms: this is a matrix exponential Fourier transform based on Ell's original two-sided 
two-dimensional quaternion transform [51 Theorem 4.1], [7J, |27j . A more general formulation is: 

Af-l JV-l 

F[u,v] = SJ2 ^e- J2v ^ f[m,n]e- Ki ^ (5) 

rn— n— 
Af-l N-l 

f[m, n] = T E e +J2 **w F[u, v}e +K27T ^ (6) 
U =Q v=a 

in which each element of the two-dimensional arrays F and / is a square matrix representing a complex or 
hypercomplex number using a matrix isomorphism for the algebra in use, for example the representations 
already given in § l4.2l in the case of the quaternion algebra; the two scale factors multiply to give 1/MTV, and 
J and K are matrix representations of two arbitrary roots of —1 in the chosen algebra. (In Ell's original 
formulation, the roots of —1 were j and k, that is two of the orthogonal quaternion basis elements. The 
following theorem shows that there is no requirement for the two roots to be orthogonal in order for the 
transform to invert.) 

Theorem 3. The transforms in ([5]) and ^ are a two-dimensional discrete Fourier transform pair, provided 
that J 2 = K 2 = I. 

Proof. The proof follows the same scheme as the proof of Theorem [TJ but we adopt a more concise presen- 
tation to fit the available column space. We start by substituting @ into , replacing m and n by M. and 
M respectively to keep the indices distinct: 

u=0 u=0 
Af-l N-l 

_A4=0Af=0 

x e K2 ^ 

The scale factors can be moved outside both summations, and replaced with their product 1/MTV; and the 
exponentials of the outer summations can be moved inside the inner, because they are constant with respect 
to the summation indices M. and M . At the same time, adjacent exponentials with the same root of — 1 
can be merged. With these changes, and omitting the scale factor to save space, the right-hand side of the 
equation becomes: 

Af-l N-l A/-1 Af-l 

E E E £ e^^/IM^e^ 2 ^ 

u=0 v=0 M=0Af=Q 

We now isolate out from the inner pair of summations the case where M = m and M = n. In this case the 
exponentials reduce to identity matrices, and we have: 

Af-liV-l 



MTV 

u=0 v=0 
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This sums to f[m, n], which is the original two-dimensional signal, as required. To complete the proof we have 
to show that the rest of the summation, excluding the case M. = m and Af = n, reduces to zero. Dropping 
the scale factor, and changing the order of summation, we have the following inner double summation: 

M-l JV-l 
u=0 v=Q 

Noting that the first exponential and / arc independent of the second summation index v, we can move 
them outside the second summation (we could do similarly with the exponential on the right and the first 
summation): 

M-l N-l 

]T e^ 1 ^^ f[M,Af] £ e K ^^^ 

and, as in Theorem [T] the summation on the right is over an integral number of cycles of cosine and sine, 
and therefore vanishes. □ □ 

Notice that it was not necessary to assume that J and K were orthogonal: it is sufficient that each be 
a root of —1. This has been verified numerically using the two-dimensional code given in the Appendix. 



8 Discussion 

We have shown that any discrete Fourier transform in an algebra that has a matrix representation, can 
be formulated in the way shown here. This includes the complex, quaternion, biquaternion, and Clifford 
algebras (although we have demonstrated only certain cases of Clifford algebras, we believe the result holds 
in general). This observation provides a theoretical unification of diverse hypercomplex DFTs. 

Several immediate possibilities for further work, as well as ramifications, now suggest themselves. Firstly, 
the study of roots of —1 is accessible from the matrix representation as well as direct representation in 
whatever algebra is employed for the transform. All of the results obtained so far in hypercomplex algebras, 
and known to the authors [25( pp203, 209], [17 ( I19j. were achieved by working in the algebra in question, 
that is by algebraic manipulation of quaternion, biquaternion or Clifford multivector values. An alternative 
approach would be to work in the equivalent matrix algebra, but this seems difficult even for the lower order 
cases. Nevertheless, it merits further study because of the possibility of finding a systematic approach that 
would cover many algebras in one framework. Following the reasoning in §[5j it is possible to define matrix 
roots of —1 that appear not to be isomorphic to any Clifford or quaternion algebra, and these merit further 
study. 

Secondly, the matrix formulation presented here lends itself to analysis of the structure of the transform, 
including possible factorizations for fast algorithms, as well as parallel or vectorized implementations for 
single-instruction, multiple-data (simd) processors, and of course, factorizations into multiple complex FFTs 
as has been done for quaternion FFTs (see for example |15j). In the case of matrix roots of —1 which do 
not correspond to Clifford or quaternion algebras, analysis of the structure of the transform may give insight 
into possible applications of transforms based on such roots. 

Finally, at a practical level, hypercomplex transforms implemented directly in hypercomplex arithmetic 
are likely to be much faster than any implementation based on matrices, but the simplicity of the matrix 
exponential formulation discussed in this paper, and the fact that it can be computed using standard real or 
complex matrix arithmetic, without using a hypercomplex library, means that the matrix exponential formu- 
lation provides a very simple reference implementation which can be used for verification of the correctness of 
hypercomplex code. This is an important point, because verification of the correctness of hypercomplex FFT 
code is otherwise non-trivial. Verification of inversion is simple enough, but establishing that the spectral 
coefficients have the correct values is much less so. 
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A mat lab® code 



We include here two short matlab® functions for computing the forward transform given in ([2]), and (|5|), 
apart from the scale factors. The inverses can be computed simply by interchanging the input and output 
and negating the matrix roots of —1. Neither function is coded for speed, on the contrary the coding is 
intended to be simple and easily verified against the equations. 

function F = matdft(f , J) 
M = size(f , 2) ; 
F = zeros (size (f) ) ; 
for m = 0:M-1 
for u = 0:M-1 

F(: , u + 1) = F(: , u + 1) ... 

+ expm(-J .* 2 .* pi .* m .* u./M) ... 

* f (: , m + 1) ; 
end 
end 

function F = matdft2(f , J, K) 
A = size (J , 1) ; 
M = size(f , 1) ./ A; 
N = size(f, 2) ./ A; 
F = zeros (size (f) ) ; 
for u = 0:M-1 
for v = 0:N-1 
for m = 0:M-1 
for n = 0:N-1 

F(A*u+l : A*u+A, A*v+l:A*v+A) = ... 
F(A*u+l : A*u+A, A*v+l:A*v+A) + ... 
expm(-J .* 2*pi .* m .* u./M) ... 

* f (A*m+1 : A*m+A, A*n+l:A*n+A) ... 

* expm(-K .* 2*pi . * n .* v./N); 
end 

end 
end 
end 
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